﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "b_bessel_1")]
    public static unsafe double b_bessel_1_dll(int n, double x)
    {
        return b_bessel_1(n, x);
    }

    /// <summary>
    /// 变型第一类整数阶Bessel函数
    /// </summary>
    /// <param name="n">n阶数。要求n>0。</param>
    /// <param name="x">x自变量值。</param>
    /// <returns>函数返回变型第一类整数阶Bessel函数值。</returns>
    public static unsafe double b_bessel_1(int n, double x)
    {
        int i, m;
        double t, y, p = 0.0, b0, b1, q;
        double[] a = new double[7] { 1.0, 3.5156229, 3.0899424, 1.2067492, 0.2659732, 0.0360768, 0.0045813 };
        double[] b = new double[7] { 0.5, 0.87890594, 0.51498869, 0.15084934, 0.02658773, 0.00301532, 0.00032411 };
        double[] c = new double[9] { 0.39894228, 0.01328592, 0.00225319, -0.00157565, 0.00916281, -0.02057706, 0.02635537, -0.01647633, 0.00392377 };
        double[] d = new double[9] { 0.39894228, -0.03988024, -0.00362018, 0.00163801, -0.01031555, 0.02282967, -0.02895312, 0.01787654, -0.00420059 };

        if (n < 0) n = -n;
        t = Math.Abs(x);
        if (n != 1)
        {
            if (t < 3.75)
            {
                y = (x / 3.75) * (x / 3.75); p = a[6];
                for (i = 5; i >= 0; i--) p = p * y + a[i];
            }
            else
            {
                y = 3.75 / t; p = c[8];
                for (i = 7; i >= 0; i--) p = p * y + c[i];
                p = p * Math.Exp(t) / Math.Sqrt(t);
            }
        }
        if (n == 0) return (p);
        q = p;
        if (t < 3.75)
        {
            y = (x / 3.75) * (x / 3.75); p = b[6];
            for (i = 5; i >= 0; i--) p = p * y + b[i];
            p = p * t;
        }
        else
        {
            y = 3.75 / t; p = d[8];
            for (i = 7; i >= 0; i--) p = p * y + d[i];
            p = p * Math.Exp(t) / Math.Sqrt(t);
        }
        if (x < 0.0) p = -p;
        if (n == 1) return (p);
        if (x == 0.0) return (0.0);
        y = 2.0 / t; t = 0.0; b1 = 1.0; b0 = 0.0;
        m = n + (int)Math.Sqrt(40.0 * n);
        m = 2 * m;
        for (i = m; i > 0; i--)
        {
            p = b0 + i * y * b1; b0 = b1; b1 = p;
            if (Math.Abs(b1) > 1.0e+10)
            {
                t = t * 1.0e-10; b0 = b0 * 1.0e-10;
                b1 = b1 * 1.0e-10;
            }
            if (i == n) t = b0;
        }
        p = t * q / b1;
        if ((x < 0.0) && (n % 2 == 1)) p = -p;
        return (p);
    }

    /*
    // 变型第一类整数阶Bessel函数例
      int main()
      { 
          int n,i;
          double x,y;
          for (n=0; n<=5; n++)
          { 
              x=0.05;
              for (i=1; i<=4; i++)
              { 
                  y=b_bessel_1(n,x);
                  cout <<"n=" <<n <<"   x=" <<setw(5) <<x 
                       <<"    I(n,x)=" <<y <<endl;
                  x=x*10.0;
              }
          }
          return 0;
      }
    */
}

